Bee Respiration for Susie and Stacey
Callin Switzer
11/4/2016, Update: 12/13/2016, Update: 28 Jan 2017
Update 28 Jan 2017
Summary of findings:
- % load and total mass are too collinear to include in the same model.
- Treatment and % load are too collinear to include in the same model.
- For the remainder of the modeling, I’ll use only % load
- Conducted all tests for squared terms of continuous predictors (suspect wbf matters)
Differences between model with Treatment vs. model with % load
Treatment order is not a significant predictor of wingbeat frequency in the model with percent load (though it is very close, p ~ 0.07)
Percent load and treatment suggest very similar models in terms of wingbeat frequency
Treatment order is not a strong predictor of stroke amplitude in the model with percent load
wingbeat freq = size Met rate = mass carrying amplitude = percent loading
Susie’s (hypotheses): lm( frequency ~ ITspan + trt order)
lm( metabolism ~ total mass)
lm( stroke amplitude ~ %load + trt) (those are probably the same thing)
- Redo the model with %load and/or total mass (probably the best predictor of metabolism)
- Think about how to make figures
How heavy the bee is loaded overall, affects how much the amplitude and frequency change
Questions to look into: ** A really heavy loaded bee will have a smaller change in metabolism vs. light bee ** Suggests a non-linear relationship or random slopes that depend on bee size — the cost for additional weight goes down as you get heavier loaded. — not the case for the change in amplitude (change in load vs. change in amplitude is linear) — not the case for metabolism or wingbeat frequency.
** Hypotheses: 1. Wingbeat freq moves around for a number of reasons (tired, cold, added mass) 2. Stroke amplitude rescues wingbeat frequency 3. As you get really close to your max, adding more weight doesn’t seem to cost as much.
Analyses of respirometry data:
- How does carrying a load during flight affect bumblebees’?
- How does load affect respiratory rate?
- How does load affect wingbeat frequency?
- How does load affect stroke amplitude?
- How does load affect wing velocity?
- Which flight measurements are most closely associated with respiratory rate in bees – wingbeat frequency, stroke amplitude, or wing velocity?
Setup
Install required packages and read in data Define custom function for evaluating VIF with multilevel models
ipak <- function(pkg){
new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
if(length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
sapply(pkg, require, character.only = TRUE)
}
packages <- c('ggplot2', 'car', 'lme4', 'gsheet', "MASS", 'influence.ME', 'sjPlot')
ipak(packages)## ggplot2 car lme4 gsheet MASS
## TRUE TRUE TRUE TRUE TRUE
## influence.ME sjPlot
## TRUE TRUEtheme_set(theme_bw())
# read in data -- google sheet called "Bumble mumble grumble"
bdta <- gsheet2tbl("https://docs.google.com/spreadsheets/d/1GUUoFq41Ep3sJNFXiyMcp7fZhYmQCzTcLnhVXdr4WRo/edit?usp=sharing")Function for calculating vif for lmer
(used later)
vif.mer <- function (fit) {
## adapted from rms::vif
v <- vcov(fit)
nam <- names(fixef(fit))
## exclude intercepts
ns <- sum(1 * (nam == "Intercept" | nam == "(Intercept)"))
if (ns > 0) {
v <- v[-(1:ns), -(1:ns), drop = FALSE]
nam <- nam[-(1:ns)]
}
d <- diag(v)^0.5
v <- diag(solve(v/(d %o% d)))
names(v) <- nam
v
}Data Overview and Engineering
# create centered and squared variables
bdta <- within(bdta, {
# new variables
percLoad <- (Mstarved + load) / Mstarved
totalMass <- Mstarved + load})
bdta <- within(bdta, {
# centered variables
percLoad_cent <- scale(bdta$percLoad, center = TRUE, scale = FALSE)[,1]
totalMass_cent <- scale(totalMass, center = TRUE, scale = FALSE)[,1]
wbf_cent <- scale(wbf.aud., center = TRUE, scale = FALSE)[,1]
strokeAmp_cent <- scale(stroke.amplitude, center = TRUE, scale = FALSE)[,1]
wingVel_cent <- scale(wing.velocity, center = TRUE, scale = FALSE)[,1]
#squared
percLoad2 <- percLoad_cent^2
totalMass2 <- totalMass^2
wbf2 <- wbf_cent^2
strokeAmp2 <- strokeAmp_cent^2
wingVel2 <- wingVel_cent^2
# change factor properties
# convert trt order to factor
Treatment.order <- as.factor(as.character(Treatment.order))
# change reference level to unloaded
Treatment <- factor(Treatment, levels = c("UL", "L"))
})
summary(bdta)## Bee.ID Treatment.order Treatment Mstarved
## Length:60 1:30 UL:30 Min. :0.0767
## Class :character 2:30 L :30 1st Qu.:0.1213
## Mode :character Median :0.1378
## Mean :0.1411
## 3rd Qu.:0.1670
## Max. :0.1909
## load IT.Span single.wing.area av.resp..CO2.mL.hr.
## Min. :0.01650 Min. :3.860 Min. :20.03 Min. : 5.507
## 1st Qu.:0.04713 1st Qu.:4.400 1st Qu.:27.93 1st Qu.: 8.935
## Median :0.07840 Median :4.620 Median :30.67 Median :10.666
## Mean :0.08332 Mean :4.594 Mean :30.70 Mean :10.600
## 3rd Qu.:0.11800 3rd Qu.:4.870 3rd Qu.:33.75 3rd Qu.:12.297
## Max. :0.17460 Max. :5.180 Max. :39.05 Max. :16.475
## wbf.aud. stroke.amplitude wing.velocity totalMass
## Min. :163.9 Min. : 96.88 Min. :1.208 Min. :0.1062
## 1st Qu.:173.0 1st Qu.:114.16 1st Qu.:1.351 1st Qu.:0.1808
## Median :179.1 Median :125.96 Median :1.427 Median :0.2191
## Mean :178.9 Mean :123.31 Mean :1.465 Mean :0.2245
## 3rd Qu.:185.3 3rd Qu.:133.61 3rd Qu.:1.588 3rd Qu.:0.2646
## Max. :194.7 Max. :144.78 Max. :1.843 Max. :0.3550
## percLoad wingVel2 strokeAmp2
## Min. :1.143 Min. :7.978e-05 Min. : 0.0093
## 1st Qu.:1.321 1st Qu.:3.980e-03 1st Qu.: 42.9577
## Median :1.576 Median :1.444e-02 Median : 99.2962
## Mean :1.600 Mean :2.509e-02 Mean :148.2959
## 3rd Qu.:1.863 3rd Qu.:3.750e-02 3rd Qu.:188.1689
## Max. :2.180 Max. :1.423e-01 Max. :698.3210
## wbf2 totalMass2 percLoad2
## Min. : 0.00127 Min. :0.01128 Min. :0.0003255
## 1st Qu.: 6.59393 1st Qu.:0.03270 1st Qu.:0.0353246
## Median : 38.46046 Median :0.04801 Median :0.0764175
## Mean : 57.11835 Mean :0.05374 Mean :0.0964447
## 3rd Qu.: 77.70129 3rd Qu.:0.07003 3rd Qu.:0.1477398
## Max. :250.38161 Max. :0.12602 Max. :0.3363085
## wingVel_cent strokeAmp_cent wbf_cent
## Min. :-0.25706 Min. :-26.426 Min. :-15.074
## 1st Qu.:-0.11484 1st Qu.: -9.145 1st Qu.: -5.972
## Median :-0.03848 Median : 2.656 Median : 0.185
## Mean : 0.00000 Mean : 0.000 Mean : 0.000
## 3rd Qu.: 0.12225 3rd Qu.: 10.307 3rd Qu.: 6.394
## Max. : 0.37723 Max. : 21.473 Max. : 15.823
## totalMass_cent percLoad_cent
## Min. :-0.118253 Min. :-0.45679
## 1st Qu.:-0.043628 1st Qu.:-0.27912
## Median :-0.005353 Median :-0.02372
## Mean : 0.000000 Mean : 0.00000
## 3rd Qu.: 0.040172 3rd Qu.: 0.26297
## Max. : 0.130547 Max. : 0.57992ggplot(bdta, aes(x = percLoad, fill = Treatment)) +
geom_histogram() ## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# seems like percent load and trt are very related# visualize % load vs total mass
ggplot(bdta, aes(x = percLoad, y = totalMass, color = Treatment)) +
geom_point()
# seem to be strongly associated# see number of observations per bee
table(bdta$Bee.ID) # each bee has two observations##
## E01 E03 E05 E09 E10 E11 E12 E16 E20 E24 E26 E28 E30 E31 E32 E33 E35 E36
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## E37 E39 E41 E42 E44 E45 E49 E50 E53 E54 E56 E59
## 2 2 2 2 2 2 2 2 2 2 2 2# get number of unque bees
length(table(bdta$Bee.ID))## [1] 30# get treatment orders
trtOrders_loadedSecond <- sapply(X = unique(bdta$Bee.ID), FUN = function(x){
tmp = bdta[bdta$Bee.ID == x, c("Treatment.order", "Treatment")]
if("2_L" %in% paste(tmp[,1], tmp[,2], sep = "_")){
loadedSecond = TRUE
}
else loadedSecond = FALSE
return(loadedSecond)
})
table(trtOrders_loadedSecond)## trtOrders_loadedSecond
## FALSE
## 30# visualize scatterplot
car::scatterplotMatrix(bdta[, 4:13])
Variance Inflation Factor and PCA
# VIF is high among the bee size predictors
# note: this function is almost the same as car::vif
vif.mer(lmer(av.resp..CO2.mL.hr. ~ Treatment + Treatment.order + Mstarved + IT.Span + totalMass + single.wing.area + percLoad + (1|Bee.ID), data = bdta))## TreatmentL Treatment.order2 Mstarved IT.Span
## 16.591645 1.270430 23.112962 6.658432
## totalMass single.wing.area percLoad
## 32.931859 10.903672 38.415618car::scatterplotMatrix(bdta[, c("Mstarved", "IT.Span", "single.wing.area")])
# principle components
aa = prcomp(bdta[, c("Mstarved", "IT.Span", "single.wing.area")], center = TRUE, scale = TRUE)
summary(aa) # 1st pc explains ~95% of the variance in the three measurements of size## Importance of components:
## PC1 PC2 PC3
## Standard deviation 1.6852 0.33843 0.21323
## Proportion of Variance 0.9467 0.03818 0.01516
## Cumulative Proportion 0.9467 0.98484 1.00000biplot(aa) # shows that all three size measurement are correlated
# note, I changed the signs of the predictions so that higher PC1 values
# correspond to bigger bees
p1 = -predict(aa)[,1]
# add PC1 scores to dataset
bdta$size_pc1 = p1
# show scatterplot matrix to see correlations among size predictors
car::scatterplotMatrix(bdta[, c("Mstarved", "IT.Span", "single.wing.area", "size_pc1", "percLoad")])
# check VIF one more time
# VIF is high among the bee size predictors
vif.mer(lmer(av.resp..CO2.mL.hr. ~ Treatment + Treatment.order + size_pc1 + percLoad + totalMass + (1|Bee.ID), data = bdta))## TreatmentL Treatment.order2 size_pc1 percLoad
## 16.643643 1.258990 7.448179 30.631628
## totalMass
## 26.790114# remove treatment
vif.mer(lmer(av.resp..CO2.mL.hr. ~ Treatment.order + size_pc1 + percLoad + totalMass + (1|Bee.ID), data = bdta))## Treatment.order2 size_pc1 percLoad totalMass
## 1.014834 7.698050 20.910198 25.813622# looks like we cannot have total mass and percentLoad in the same model, according to VIF
# remove total mass
vif.mer(lmer(av.resp..CO2.mL.hr. ~ Treatment.order + size_pc1 + percLoad + (1|Bee.ID), data = bdta))## Treatment.order2 size_pc1 percLoad
## 1.003811 1.004814 1.008624Multilevel models
Resp Rate
# make a full model with all two-way interactions
m1 <- lmer(av.resp..CO2.mL.hr. ~ (size_pc1 + Treatment.order + percLoad_cent)^2 +
percLoad2 +
+ (1|Bee.ID), data = bdta)
summary(m1)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## av.resp..CO2.mL.hr. ~ (size_pc1 + Treatment.order + percLoad_cent)^2 +
## percLoad2 + +(1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 182.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7335 -0.4895 0.1150 0.4689 2.3779
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.4520 1.2050
## Residual 0.5021 0.7086
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 11.23228 0.27961 40.17
## size_pc1 1.28130 0.15748 8.14
## Treatment.order2 -0.28829 0.18935 -1.52
## percLoad_cent 4.90035 0.70897 6.91
## percLoad2 -5.10731 1.38346 -3.69
## size_pc1:Treatment.order2 -0.17569 0.12108 -1.45
## size_pc1:percLoad_cent 0.09533 0.21461 0.44
## Treatment.order2:percLoad_cent -1.19963 1.15631 -1.04
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2 s_1:T. s_1:L_
## size_pc1 -0.001
## Trtmnt.rdr2 -0.253 0.032
## percLod_cnt -0.042 0.233 0.066
## percLoad2 -0.392 -0.049 -0.204 -0.043
## sz_pc1:Tr.2 -0.049 -0.410 -0.089 -0.325 0.219
## sz_pc1:prL_ -0.030 -0.105 -0.205 -0.146 0.331 0.289
## Trtmnt.2:L_ 0.071 -0.207 -0.029 -0.886 -0.035 0.284 0.159### LRT's for interactions
m2.0 <- update(m1, .~. - size_pc1:percLoad_cent)
anova(m1, m2.0) # likelihood ratio test for interaction of treatment order## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m2.0: av.resp..CO2.mL.hr. ~ size_pc1 + Treatment.order + percLoad_cent +
## m2.0: percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order + Treatment.order:percLoad_cent
## m1: av.resp..CO2.mL.hr. ~ (size_pc1 + Treatment.order + percLoad_cent)^2 +
## m1: percLoad2 + +(1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.0 9 195.92 214.77 -88.961 177.92
## m1 10 197.68 218.62 -88.841 177.68 0.2397 1 0.6244summary(m2.0)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## av.resp..CO2.mL.hr. ~ size_pc1 + Treatment.order + percLoad_cent +
## percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order + Treatment.order:percLoad_cent
## Data: bdta
##
## REML criterion at convergence: 181.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.72009 -0.52773 0.07283 0.48826 2.35677
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.4603 1.2084
## Residual 0.4862 0.6973
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 11.2371 0.2783 40.37
## size_pc1 1.2881 0.1562 8.24
## Treatment.order2 -0.2709 0.1824 -1.49
## percLoad_cent 4.9339 0.6939 7.11
## percLoad2 -5.3211 1.2854 -4.14
## size_pc1:Treatment.order2 -0.1907 0.1141 -1.67
## Treatment.order2:percLoad_cent -1.2596 1.1316 -1.11
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2 s_1:T.
## size_pc1 -0.004
## Trtmnt.rdr2 -0.262 0.010
## percLod_cnt -0.048 0.219 0.036
## percLoad2 -0.401 -0.014 -0.147 0.008
## sz_pc1:Tr.2 -0.041 -0.393 -0.031 -0.299 0.136
## Trtmnt.2:L_ 0.078 -0.192 0.004 -0.885 -0.096 0.253m2.1 <- update(m2.0,.~. - Treatment.order:percLoad_cent)
anova(m2.0, m2.1)## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m2.1: av.resp..CO2.mL.hr. ~ size_pc1 + Treatment.order + percLoad_cent +
## m2.1: percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order
## m2.0: av.resp..CO2.mL.hr. ~ size_pc1 + Treatment.order + percLoad_cent +
## m2.0: percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order + Treatment.order:percLoad_cent
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.1 8 195.15 211.90 -89.573 179.15
## m2.0 9 195.92 214.77 -88.961 177.92 1.2255 1 0.2683summary(m2.1)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## av.resp..CO2.mL.hr. ~ size_pc1 + Treatment.order + percLoad_cent +
## percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order
## Data: bdta
##
## REML criterion at convergence: 185
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.69131 -0.46047 0.00547 0.49254 2.34469
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.5578 1.2481
## Residual 0.4642 0.6813
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 11.2638 0.2811 40.07
## size_pc1 1.2546 0.1562 8.03
## Treatment.order2 -0.2696 0.1782 -1.51
## percLoad_cent 4.2487 0.3162 13.44
## percLoad2 -5.4887 1.2518 -4.38
## size_pc1:Treatment.order2 -0.1589 0.1079 -1.47
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2
## size_pc1 0.010
## Trtmnt.rdr2 -0.254 0.011
## percLod_cnt 0.044 0.102 0.086
## percLoad2 -0.383 -0.032 -0.147 -0.165
## sz_pc1:Tr.2 -0.061 -0.348 -0.033 -0.167 0.167m2.2 <- update(m2.1, .~. - size_pc1:Treatment.order)
anova(m2.1, m2.2) ## drop all interactions## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m2.2: av.resp..CO2.mL.hr. ~ size_pc1 + Treatment.order + percLoad_cent +
## m2.2: percLoad2 + (1 | Bee.ID)
## m2.1: av.resp..CO2.mL.hr. ~ size_pc1 + Treatment.order + percLoad_cent +
## m2.1: percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.2 7 195.54 210.21 -90.773 181.54
## m2.1 8 195.15 211.90 -89.573 179.15 2.3985 1 0.1215summary(m2.2)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## av.resp..CO2.mL.hr. ~ size_pc1 + Treatment.order + percLoad_cent +
## percLoad2 + (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 184.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6892 -0.4659 0.0722 0.3830 2.1677
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.5420 1.2418
## Residual 0.4856 0.6968
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 11.2373 0.2818 39.88
## size_pc1 1.1746 0.1462 8.03
## Treatment.order2 -0.2786 0.1822 -1.53
## percLoad_cent 4.1717 0.3187 13.09
## percLoad2 -5.1672 1.2616 -4.10
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 prcLd_
## size_pc1 -0.012
## Trtmnt.rdr2 -0.261 -0.001
## percLod_cnt 0.035 0.049 0.081
## percLoad2 -0.385 0.029 -0.144 -0.141# renaming model to simplify later typing
m2 <- m2.2
##### LRTs for main effects
## Treatment Order
m3 <- update(m2, .~. - Treatment.order)
anova(m2, m3, test = "Chi") # drop trt order (different than model with treatment)## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m3: av.resp..CO2.mL.hr. ~ size_pc1 + percLoad_cent + percLoad2 +
## m3: (1 | Bee.ID)
## m2: av.resp..CO2.mL.hr. ~ size_pc1 + Treatment.order + percLoad_cent +
## m2: percLoad2 + (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m3 6 196.04 208.60 -92.019 184.04
## m2 7 195.54 210.21 -90.773 181.54 2.4925 1 0.1144summary(m3)## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ size_pc1 + percLoad_cent + percLoad2 +
## (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 185.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.61825 -0.42513 -0.04217 0.52054 2.29988
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.5364 1.2395
## Residual 0.5073 0.7123
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 11.1235 0.2735 40.67
## size_pc1 1.1745 0.1465 8.02
## percLoad_cent 4.2122 0.3245 12.98
## percLoad2 -5.4323 1.2754 -4.26
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 prcLd_
## size_pc1 -0.013
## percLod_cnt 0.059 0.050
## percLoad2 -0.450 0.029 -0.131# LRT for size
m4 <- update(m3, .~. - percLoad2)
anova(m3, m4) # keep percLoad2## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m4: av.resp..CO2.mL.hr. ~ size_pc1 + percLoad_cent + (1 | Bee.ID)
## m3: av.resp..CO2.mL.hr. ~ size_pc1 + percLoad_cent + percLoad2 +
## m3: (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m4 5 208.25 218.73 -99.127 198.25
## m3 6 196.04 208.60 -92.019 184.04 14.216 1 0.0001629 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m3)## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ size_pc1 + percLoad_cent + percLoad2 +
## (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 185.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.61825 -0.42513 -0.04217 0.52054 2.29988
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.5364 1.2395
## Residual 0.5073 0.7123
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 11.1235 0.2735 40.67
## size_pc1 1.1745 0.1465 8.02
## percLoad_cent 4.2122 0.3245 12.98
## percLoad2 -5.4323 1.2754 -4.26
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 prcLd_
## size_pc1 -0.013
## percLod_cnt 0.059 0.050
## percLoad2 -0.450 0.029 -0.131# LRT for Treatment (load)
m5 <- update(m3, .~. - size_pc1)
anova(m3, m5) # keep size## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m5: av.resp..CO2.mL.hr. ~ percLoad_cent + percLoad2 + (1 | Bee.ID)
## m3: av.resp..CO2.mL.hr. ~ size_pc1 + percLoad_cent + percLoad2 +
## m3: (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m5 5 229.84 240.31 -109.920 219.84
## m3 6 196.04 208.60 -92.019 184.04 35.801 1 2.185e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m6 <- update(m3, .~. - percLoad_cent)
anova(m3, m6) # keep squared term## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m6: av.resp..CO2.mL.hr. ~ size_pc1 + percLoad2 + (1 | Bee.ID)
## m3: av.resp..CO2.mL.hr. ~ size_pc1 + percLoad_cent + percLoad2 +
## m3: (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m6 5 257.02 267.49 -123.510 247.02
## m3 6 196.04 208.60 -92.019 184.04 62.982 1 2.086e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# summarize final model for paper
summary(m3)## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ size_pc1 + percLoad_cent + percLoad2 +
## (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 185.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.61825 -0.42513 -0.04217 0.52054 2.29988
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.5364 1.2395
## Residual 0.5073 0.7123
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 11.1235 0.2735 40.67
## size_pc1 1.1745 0.1465 8.02
## percLoad_cent 4.2122 0.3245 12.98
## percLoad2 -5.4323 1.2754 -4.26
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 prcLd_
## size_pc1 -0.013
## percLod_cnt 0.059 0.050
## percLoad2 -0.450 0.029 -0.131# write output
summary(m3)$coefficients ## Estimate Std. Error t value
## (Intercept) 11.123542 0.2734947 40.671875
## size_pc1 1.174509 0.1464512 8.019799
## percLoad_cent 4.212203 0.3245350 12.979195
## percLoad2 -5.432315 1.2754342 -4.259189write.csv( summary(m3)$coefficients, file = "RespCoefs_percLoad.csv" )resp rate diagnostics
# qq plot
qqnorm(resid(m3), main = "")
qqline(resid(m3)) # good
# residual plot
plot(fitted(m3), resid(m3), xlab = "fitted", ylab = "residuals")
abline(0,0)
# check cook's distance -- should be less than 1, so we're good
infl <- influence(m3, obs = TRUE)
plot(infl, which = 'cook')
# QQPlot for group-level effects
qqnorm(ranef(m3)$Bee.ID[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(m3)$Bee.ID[[1]]) # looks good
resp rate visualization
newdf <- data.frame(size_pc1 = 0, Treatment.order = factor(1), percLoad_cent = seq(min(bdta$percLoad_cent), max(bdta$percLoad_cent), length = 100))
newdf$percLoad2 <- newdf$percLoad_cent^2
preds1 <- predict(m3, newdf, re.form = NA)
# plot of prediction for average sized bee (line)
# raw data plotted as points
ggplot(bdta, aes(x = percLoad_cent + mean(bdta$percLoad) , y = av.resp..CO2.mL.hr.)) +
geom_point(alpha = 0.3) +
geom_line(data = newdf, aes(x = percLoad_cent + mean(bdta$percLoad), y = preds1)) +
labs(x = "Mass during flight (% of starved mass in g)")
Wingbeat Freq
# make a full model with all two-way interactions
m1 <- lmer(wbf.aud. ~ (size_pc1 + Treatment.order + percLoad_cent)^2 +
percLoad2 +
(1|Bee.ID), data = bdta)
summary(m1)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## wbf.aud. ~ (size_pc1 + Treatment.order + percLoad_cent)^2 + percLoad2 +
## (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 326.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8304 -0.7064 0.1416 0.5616 1.5135
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 17.295 4.159
## Residual 9.876 3.143
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 183.2406 1.0745 170.53
## size_pc1 -2.9657 0.5916 -5.01
## Treatment.order2 -2.6627 0.8394 -3.17
## percLoad_cent 7.0660 2.9186 2.42
## percLoad2 -31.4157 6.0747 -5.17
## size_pc1:Treatment.order2 0.2089 0.5339 0.39
## size_pc1:percLoad_cent -0.5715 0.9461 -0.60
## Treatment.order2:percLoad_cent -0.7048 4.6314 -0.15
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2 s_1:T. s_1:L_
## size_pc1 0.005
## Trtmnt.rdr2 -0.294 0.038
## percLod_cnt -0.026 0.256 0.072
## percLoad2 -0.446 -0.066 -0.202 -0.081
## sz_pc1:Tr.2 -0.062 -0.476 -0.089 -0.311 0.228
## sz_pc1:prL_ -0.035 -0.121 -0.204 -0.146 0.332 0.286
## Trtmnt.2:L_ 0.061 -0.224 -0.034 -0.868 0.001 0.267 0.161### LRT's for interactions
m2.0 <- update(m1, .~. - Treatment.order:percLoad_cent)
anova(m1, m2.0)## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m2.0: wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 +
## m2.0: (1 | Bee.ID) + size_pc1:Treatment.order + size_pc1:percLoad_cent
## m1: wbf.aud. ~ (size_pc1 + Treatment.order + percLoad_cent)^2 + percLoad2 +
## m1: (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.0 9 362.18 381.03 -172.09 344.18
## m1 10 364.15 385.10 -172.08 344.15 0.0311 1 0.86summary(m2.0)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 +
## (1 | Bee.ID) + size_pc1:Treatment.order + size_pc1:percLoad_cent
## Data: bdta
##
## REML criterion at convergence: 331.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8513 -0.6994 0.1341 0.5629 1.5305
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 16.921 4.113
## Residual 9.716 3.117
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 183.2499 1.0623 172.50
## size_pc1 -2.9860 0.5709 -5.23
## Treatment.order2 -2.6673 0.8321 -3.21
## percLoad_cent 6.6792 1.4349 4.65
## percLoad2 -31.4059 6.0246 -5.21
## size_pc1:Treatment.order2 0.2308 0.5103 0.45
## size_pc1:percLoad_cent -0.5481 0.9261 -0.59
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2 s_1:T.
## size_pc1 0.020
## Trtmnt.rdr2 -0.293 0.031
## percLod_cnt 0.054 0.128 0.086
## percLoad2 -0.447 -0.068 -0.202 -0.162
## sz_pc1:Tr.2 -0.082 -0.444 -0.082 -0.164 0.237
## sz_pc1:prL_ -0.045 -0.088 -0.201 -0.012 0.336 0.255m2.1 <- update(m2.0,.~. - size_pc1:Treatment.order)
anova(m2.0, m2.1)## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m2.1: wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 +
## m2.1: (1 | Bee.ID) + size_pc1:percLoad_cent
## m2.0: wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 +
## m2.0: (1 | Bee.ID) + size_pc1:Treatment.order + size_pc1:percLoad_cent
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.1 8 360.42 377.17 -172.21 344.42
## m2.0 9 362.18 381.03 -172.09 344.18 0.2349 1 0.6279summary(m2.1)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 +
## (1 | Bee.ID) + size_pc1:percLoad_cent
## Data: bdta
##
## REML criterion at convergence: 332.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8110 -0.6782 0.1540 0.5820 1.4972
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 17.141 4.140
## Residual 9.392 3.065
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 183.2949 1.0536 173.96
## size_pc1 -2.8715 0.5122 -5.61
## Treatment.order2 -2.6345 0.8154 -3.23
## percLoad_cent 6.7941 1.3926 4.88
## percLoad2 -32.1204 5.7604 -5.58
## size_pc1:percLoad_cent -0.6559 0.8809 -0.74
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2
## size_pc1 -0.018
## Trtmnt.rdr2 -0.298 -0.006
## percLod_cnt 0.041 0.062 0.074
## percLoad2 -0.437 0.042 -0.189 -0.128
## sz_pc1:prL_ -0.025 0.028 -0.187 0.031 0.294m2.2 <- update(m2.1, .~. - size_pc1:percLoad_cent)
anova(m2.1, m2.2) ## drop all interactions## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m2.2: wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 +
## m2.2: (1 | Bee.ID)
## m2.1: wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 +
## m2.1: (1 | Bee.ID) + size_pc1:percLoad_cent
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.2 7 359.05 373.71 -172.53 345.05
## m2.1 8 360.42 377.17 -172.21 344.42 0.6319 1 0.4267summary(m2.2)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 +
## (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 334.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8541 -0.7036 0.1281 0.5915 1.3890
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 17.163 4.143
## Residual 9.258 3.043
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 183.2771 1.0500 174.54
## size_pc1 -2.8608 0.5114 -5.59
## Treatment.order2 -2.7473 0.7953 -3.45
## percLoad_cent 6.8290 1.3823 4.94
## percLoad2 -30.8797 5.4673 -5.65
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 prcLd_
## size_pc1 -0.017
## Trtmnt.rdr2 -0.307 -0.001
## percLod_cnt 0.041 0.061 0.081
## percLoad2 -0.448 0.035 -0.143 -0.144# renaming model to simplify later typing
m2 <- m2.2
##### LRTs for main effects
## Treatment Order
m3 <- update(m2, .~. - Treatment.order)
anova(m2, m3, test = "Chi") # keep treatment order## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m3: wbf.aud. ~ size_pc1 + percLoad_cent + percLoad2 + (1 | Bee.ID)
## m2: wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 +
## m2: (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m3 6 368.10 380.67 -178.05 356.10
## m2 7 359.05 373.71 -172.53 345.05 11.05 1 0.0008867 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# LRT for size
m4 <- update(m2, .~. - size_pc1)
anova(m2, m4) # keep size_pc1## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m4: wbf.aud. ~ Treatment.order + percLoad_cent + percLoad2 + (1 |
## m4: Bee.ID)
## m2: wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 +
## m2: (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m4 6 379.62 392.19 -183.81 367.62
## m2 7 359.05 373.71 -172.53 345.05 22.571 1 2.025e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# LRT for Treatment (load)
m5 <- update(m2, .~. - percLoad_cent)
anova(m2, m5) # keep percLoad## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m5: wbf.aud. ~ size_pc1 + Treatment.order + percLoad2 + (1 | Bee.ID)
## m2: wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 +
## m2: (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m5 6 375.91 388.47 -181.95 363.91
## m2 7 359.05 373.71 -172.53 345.05 18.858 1 1.408e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m6 <- update(m2, .~. - percLoad2)
anova(m2, m6) # keep squared term## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m6: wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + (1 |
## m6: Bee.ID)
## m2: wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 +
## m2: (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m6 6 378.98 391.54 -183.49 366.98
## m2 7 359.05 373.71 -172.53 345.05 21.927 1 2.832e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# summarize final model for paper
summary(m2)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## wbf.aud. ~ size_pc1 + Treatment.order + percLoad_cent + percLoad2 +
## (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 334.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8541 -0.7036 0.1281 0.5915 1.3890
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 17.163 4.143
## Residual 9.258 3.043
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 183.2771 1.0500 174.54
## size_pc1 -2.8608 0.5114 -5.59
## Treatment.order2 -2.7473 0.7953 -3.45
## percLoad_cent 6.8290 1.3823 4.94
## percLoad2 -30.8797 5.4673 -5.65
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 prcLd_
## size_pc1 -0.017
## Trtmnt.rdr2 -0.307 -0.001
## percLod_cnt 0.041 0.061 0.081
## percLoad2 -0.448 0.035 -0.143 -0.144# write output
summary(m2)$coefficients ## Estimate Std. Error t value
## (Intercept) 183.277063 1.0500400 174.542933
## size_pc1 -2.860783 0.5114480 -5.593497
## Treatment.order2 -2.747317 0.7952937 -3.454469
## percLoad_cent 6.829045 1.3823068 4.940325
## percLoad2 -30.879690 5.4672717 -5.648099write.csv( summary(m2)$coefficients, file = "FreqCoefs_percLoad.csv" )wbf diagnostics
# qq plot
qqnorm(resid(m2), main = "")
qqline(resid(m2)) # good
# check influence
plot(influence(m2, obs = TRUE ), which = 'cook')
# residual plot
plot(fitted(m2), residuals(m2, "deviance"), xlab = "fitted", ylab = "residuals")
abline(0,0)
# QQPlot for group-level effects
qqnorm(ranef(m2)$Bee.ID[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(m2)$Bee.ID[[1]]) # looks good
wbf visualization
newdf <- data.frame(size_pc1 = 0, Treatment.order = factor(1), percLoad_cent = seq(min(bdta$percLoad_cent), max(bdta$percLoad_cent), length = 100))
newdf$percLoad2 <- newdf$percLoad_cent^2
preds1 <- predict(m2, newdf, re.form = NA) #predict for mean of random effect
# plot of prediction for average sized bee (line)
# raw data plotted as points
ggplot(bdta, aes(x = percLoad_cent + mean(bdta$percLoad) , y = wbf.aud.)) +
geom_point(alpha = 0.3) +
geom_line(data = newdf, aes(x = percLoad_cent + mean(bdta$percLoad), y = preds1)) +
labs(x = "Mass during flight (% of starved mass in g)")
Stroke Amplitude
# make a full model with all two-way interactions
m1 <- lmer(stroke.amplitude ~ (size_pc1 + Treatment.order + percLoad_cent)^2 +
percLoad2 +
(1|Bee.ID), data = bdta)
summary(m1)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## stroke.amplitude ~ (size_pc1 + Treatment.order + percLoad_cent)^2 +
## percLoad2 + (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 352.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.09359 -0.52135 0.01418 0.41271 2.08441
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 13.01 3.606
## Residual 24.48 4.948
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 124.1538 1.3504 91.94
## size_pc1 1.9451 0.7042 2.76
## Treatment.order2 -0.4425 1.3197 -0.34
## percLoad_cent 39.0584 3.8975 10.02
## percLoad2 -9.2460 9.3003 -0.99
## size_pc1:Treatment.order2 -1.3443 0.8311 -1.62
## size_pc1:percLoad_cent -3.1048 1.4603 -2.13
## Treatment.order2:percLoad_cent -5.6584 5.7085 -0.99
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2 s_1:T. s_1:L_
## size_pc1 0.025
## Trtmnt.rdr2 -0.374 0.049
## percLod_cnt 0.026 0.292 0.087
## percLoad2 -0.540 -0.105 -0.198 -0.167
## sz_pc1:Tr.2 -0.088 -0.611 -0.087 -0.284 0.242
## sz_pc1:prL_ -0.040 -0.151 -0.199 -0.147 0.329 0.278
## Trtmnt.2:L_ 0.016 -0.244 -0.047 -0.821 0.086 0.237 0.172### LRT's for interactions
m2.0 <- update(m1, .~. - Treatment.order:percLoad_cent )
anova(m1, m2.0)## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m2.0: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent +
## m2.0: percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order + size_pc1:percLoad_cent
## m1: stroke.amplitude ~ (size_pc1 + Treatment.order + percLoad_cent)^2 +
## m1: percLoad2 + (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.0 9 394.21 413.05 -188.10 376.21
## m1 10 395.22 416.16 -187.61 375.22 0.9852 1 0.3209summary(m2.0)## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent +
## percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order + size_pc1:percLoad_cent
## Data: bdta
##
## REML criterion at convergence: 358.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.06103 -0.56098 -0.01219 0.50169 2.15159
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 14.18 3.765
## Residual 23.65 4.864
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 124.1660 1.3486 92.07
## size_pc1 1.7778 0.6856 2.59
## Treatment.order2 -0.5058 1.2961 -0.39
## percLoad_cent 35.9655 2.1941 16.39
## percLoad2 -8.3350 9.1404 -0.91
## size_pc1:Treatment.order2 -1.1503 0.7940 -1.45
## size_pc1:percLoad_cent -2.8411 1.4174 -2.00
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2 s_1:T.
## size_pc1 0.029
## Trtmnt.rdr2 -0.367 0.038
## percLod_cnt 0.068 0.163 0.085
## percLoad2 -0.536 -0.086 -0.196 -0.168
## sz_pc1:Tr.2 -0.094 -0.575 -0.078 -0.162 0.229
## sz_pc1:prL_ -0.044 -0.112 -0.195 -0.011 0.322 0.248m2.1 <- update(m2.0,.~. - size_pc1:Treatment.order)
anova(m2.0, m2.1)## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m2.1: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent +
## m2.1: percLoad2 + (1 | Bee.ID) + size_pc1:percLoad_cent
## m2.0: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent +
## m2.0: percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order + size_pc1:percLoad_cent
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.1 8 394.61 411.36 -189.3 378.61
## m2.0 9 394.21 413.05 -188.1 376.21 2.4002 1 0.1213summary(m2.1)## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent +
## percLoad2 + (1 | Bee.ID) + size_pc1:percLoad_cent
## Data: bdta
##
## REML criterion at convergence: 362.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.02623 -0.51313 -0.04168 0.51460 2.13594
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 13.6 3.688
## Residual 24.7 4.970
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 123.9899 1.3564 91.41
## size_pc1 1.2045 0.5603 2.15
## Treatment.order2 -0.6514 1.3201 -0.49
## percLoad_cent 35.3986 2.2086 16.03
## percLoad2 -5.3965 9.0698 -0.59
## size_pc1:percLoad_cent -2.3436 1.4009 -1.67
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2
## size_pc1 -0.031
## Trtmnt.rdr2 -0.382 -0.008
## percLod_cnt 0.055 0.089 0.073
## percLoad2 -0.536 0.059 -0.183 -0.137
## sz_pc1:prL_ -0.021 0.040 -0.181 0.031 0.280m2.2 <- update(m2.1, .~. - size_pc1:percLoad_cent)
anova(m2.1, m2.2) ## drop all interactions## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m2.2: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent +
## m2.2: percLoad2 + (1 | Bee.ID)
## m2.1: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent +
## m2.1: percLoad2 + (1 | Bee.ID) + size_pc1:percLoad_cent
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.2 7 395.66 410.32 -190.83 381.66
## m2.1 8 394.61 411.36 -189.30 378.61 3.0552 1 0.08048 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m2.2)## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent +
## percLoad2 + (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 367.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.89123 -0.52057 -0.02436 0.57562 2.02535
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 13.86 3.723
## Residual 25.63 5.062
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 123.9447 1.3781 89.94
## size_pc1 1.2414 0.5677 2.19
## Treatment.order2 -1.0513 1.3225 -0.79
## percLoad_cent 35.5007 2.2478 15.79
## percLoad2 -1.1743 8.8665 -0.13
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 prcLd_
## size_pc1 -0.031
## Trtmnt.rdr2 -0.393 -0.001
## percLod_cnt 0.056 0.088 0.080
## percLoad2 -0.553 0.051 -0.140 -0.152# renaming model to simplify later typing
m2 <- m2.2
##### LRTs for main effects
## size
m3 <- update(m2, .~. - percLoad2)
anova(m2, m3, test = "Chi") # drop squared term## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m3: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent +
## m3: (1 | Bee.ID)
## m2: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent +
## m2: percLoad2 + (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m3 6 393.68 406.24 -190.84 381.68
## m2 7 395.66 410.32 -190.83 381.66 0.0174 1 0.8949summary(m3)## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent +
## (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 373.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.87059 -0.50457 -0.03627 0.58077 2.02395
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 14.31 3.782
## Residual 24.73 4.973
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 123.8429 1.1414 108.50
## size_pc1 1.2463 0.5668 2.20
## Treatment.order2 -1.0742 1.2864 -0.84
## percLoad_cent 35.5004 2.1858 16.24
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2
## size_pc1 -0.003
## Trtmnt.rdr2 -0.564 0.006
## percLod_cnt -0.034 0.096 0.060# trt order
m4 <- update(m3, .~. - Treatment.order)
anova(m3, m4) # drop treatment order## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m4: stroke.amplitude ~ size_pc1 + percLoad_cent + (1 | Bee.ID)
## m3: stroke.amplitude ~ size_pc1 + Treatment.order + percLoad_cent +
## m3: (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m4 5 392.42 402.89 -191.21 382.42
## m3 6 393.68 406.24 -190.84 381.68 0.7374 1 0.3905summary(m4)## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ size_pc1 + percLoad_cent + (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 376.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.96941 -0.48117 -0.05909 0.48300 1.94459
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 14.53 3.811
## Residual 24.42 4.941
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 123.3058 0.9440 130.62
## size_pc1 1.2495 0.5675 2.20
## percLoad_cent 35.6293 2.1693 16.42
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1
## size_pc1 0.000
## percLod_cnt 0.000 0.095# LRT for size
m5 <- update(m4, .~. - size_pc1)
anova(m4, m5) # keep size## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m5: stroke.amplitude ~ percLoad_cent + (1 | Bee.ID)
## m4: stroke.amplitude ~ size_pc1 + percLoad_cent + (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m5 4 395.24 403.62 -193.62 387.24
## m4 5 392.42 402.89 -191.21 382.42 4.828 1 0.028 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m6 <- update(m4, .~. - percLoad_cent)
anova(m4, m6) # keep percent load## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m6: stroke.amplitude ~ size_pc1 + (1 | Bee.ID)
## m4: stroke.amplitude ~ size_pc1 + percLoad_cent + (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m6 4 478.07 486.45 -235.04 470.07
## m4 5 392.42 402.89 -191.21 382.42 87.657 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# summarize final model for paper
summary(m4)## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ size_pc1 + percLoad_cent + (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 376.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.96941 -0.48117 -0.05909 0.48300 1.94459
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 14.53 3.811
## Residual 24.42 4.941
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 123.3058 0.9440 130.62
## size_pc1 1.2495 0.5675 2.20
## percLoad_cent 35.6293 2.1693 16.42
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1
## size_pc1 0.000
## percLod_cnt 0.000 0.095# write output
summary(m4)$coefficients ## Estimate Std. Error t value
## (Intercept) 123.305797 0.9440279 130.616690
## size_pc1 1.249463 0.5674547 2.201872
## percLoad_cent 35.629329 2.1693080 16.424283write.csv( summary(m4)$coefficients, file = "AmpCoefs_percLoad.csv" )strokeAmplitude diagnostics
# qq plot
qqnorm(resid(m4), main = "")
qqline(resid(m4)) # ok
# residual plot
plot(fitted(m4), resid(m3), xlab = "fitted", ylab = "residuals")
abline(0,0)
# QQPlot for group-level effects
qqnorm(ranef(m4)$Bee.ID[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(m4)$Bee.ID[[1]]) # looks good
Stroke Amplitude Visualization
preds1 <- predict(m4, newdf, re.form = NA) #predict for mean of random effect
# plot of prediction for average sized bee (line)
# raw data plotted as points
ggplot(bdta, aes(x = percLoad_cent + mean(bdta$percLoad) , y = stroke.amplitude)) +
geom_point(alpha = 0.3) +
geom_line(data = newdf, aes(x = percLoad_cent + mean(bdta$percLoad), y = preds1)) +
labs(x = "Mass during flight (% of starved mass in g)")
Wing Velocity
# make a full model with all two-way interactions
m1 <- lmer(wing.velocity ~ (size_pc1 + Treatment.order + percLoad_cent)^2 +
percLoad2 +
(1|Bee.ID), data = bdta)
summary(m1)## Linear mixed model fit by REML ['lmerMod']
## Formula: wing.velocity ~ (size_pc1 + Treatment.order + percLoad_cent)^2 +
## percLoad2 + (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: -76.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.96215 -0.44762 -0.04005 0.50444 2.11529
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 0.003968 0.06299
## Residual 0.006025 0.07762
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.48344 0.02182 67.98
## size_pc1 -0.04704 0.01148 -4.10
## Treatment.order2 -0.01763 0.02071 -0.85
## percLoad_cent -0.41560 0.06273 -6.62
## percLoad2 -0.05676 0.14666 -0.39
## size_pc1:Treatment.order2 0.01481 0.01306 1.13
## size_pc1:percLoad_cent 0.04749 0.02300 2.06
## Treatment.order2:percLoad_cent 0.05562 0.09323 0.60
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2 s_1:T. s_1:L_
## size_pc1 0.022
## Trtmnt.rdr2 -0.362 0.048
## percLod_cnt 0.017 0.288 0.085
## percLoad2 -0.527 -0.099 -0.199 -0.153
## sz_pc1:Tr.2 -0.085 -0.591 -0.087 -0.287 0.240
## sz_pc1:prL_ -0.040 -0.147 -0.200 -0.147 0.330 0.279
## Trtmnt.2:L_ 0.025 -0.242 -0.045 -0.829 0.072 0.241 0.170### LRT's for interactions
m2.0 <- update(m1, .~. - Treatment.order:percLoad_cent )
anova(m1, m2.0)## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m2.0: wing.velocity ~ size_pc1 + Treatment.order + percLoad_cent +
## m2.0: percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order + size_pc1:percLoad_cent
## m1: wing.velocity ~ (size_pc1 + Treatment.order + percLoad_cent)^2 +
## m1: percLoad2 + (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.0 9 -101.513 -82.664 59.756 -119.51
## m1 10 -99.871 -78.927 59.935 -119.87 0.3579 1 0.5497summary(m2.0)## Linear mixed model fit by REML ['lmerMod']
## Formula: wing.velocity ~ size_pc1 + Treatment.order + percLoad_cent +
## percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order + size_pc1:percLoad_cent
## Data: bdta
##
## REML criterion at convergence: -79.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.02111 -0.46259 -0.01217 0.49529 2.12092
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 0.004027 0.06346
## Residual 0.005882 0.07670
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.48322 0.02168 68.42
## size_pc1 -0.04539 0.01109 -4.09
## Treatment.order2 -0.01704 0.02044 -0.83
## percLoad_cent -0.38478 0.03470 -11.09
## percLoad2 -0.06438 0.14470 -0.44
## size_pc1:Treatment.order2 0.01292 0.01253 1.03
## size_pc1:percLoad_cent 0.04506 0.02241 2.01
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2 s_1:T.
## size_pc1 0.028
## Trtmnt.rdr2 -0.359 0.038
## percLod_cnt 0.067 0.160 0.085
## percLoad2 -0.528 -0.084 -0.197 -0.167
## sz_pc1:Tr.2 -0.093 -0.561 -0.079 -0.162 0.230
## sz_pc1:prL_ -0.045 -0.110 -0.196 -0.011 0.324 0.249m2.1 <- update(m2.0,.~. - size_pc1:Treatment.order)
anova(m2.0, m2.1)## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m2.1: wing.velocity ~ size_pc1 + Treatment.order + percLoad_cent +
## m2.1: percLoad2 + (1 | Bee.ID) + size_pc1:percLoad_cent
## m2.0: wing.velocity ~ size_pc1 + Treatment.order + percLoad_cent +
## m2.0: percLoad2 + (1 | Bee.ID) + size_pc1:Treatment.order + size_pc1:percLoad_cent
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.1 8 -102.28 -85.525 59.140 -118.28
## m2.0 9 -101.51 -82.664 59.756 -119.51 1.2331 1 0.2668summary(m2.1)## Linear mixed model fit by REML ['lmerMod']
## Formula: wing.velocity ~ size_pc1 + Treatment.order + percLoad_cent +
## percLoad2 + (1 | Bee.ID) + size_pc1:percLoad_cent
## Data: bdta
##
## REML criterion at convergence: -84.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.01405 -0.48637 0.04786 0.57874 2.09154
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 0.004062 0.06374
## Residual 0.005871 0.07662
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.485332 0.021598 68.77
## size_pc1 -0.038972 0.009194 -4.24
## Treatment.order2 -0.015368 0.020360 -0.75
## percLoad_cent -0.379047 0.034213 -11.08
## percLoad2 -0.099134 0.140706 -0.70
## size_pc1:percLoad_cent 0.039276 0.021687 1.81
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 prcLd_ prcLd2
## size_pc1 -0.029
## Trtmnt.rdr2 -0.368 -0.008
## percLod_cnt 0.052 0.084 0.073
## percLoad2 -0.522 0.056 -0.184 -0.135
## sz_pc1:prL_ -0.022 0.038 -0.182 0.031 0.283m2.2 <- update(m2.1, .~. - size_pc1:percLoad_cent)
anova(m2.1, m2.2) ## drop all interactions, though this is very close## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m2.2: wing.velocity ~ size_pc1 + Treatment.order + percLoad_cent +
## m2.2: percLoad2 + (1 | Bee.ID)
## m2.1: wing.velocity ~ size_pc1 + Treatment.order + percLoad_cent +
## m2.1: percLoad2 + (1 | Bee.ID) + size_pc1:percLoad_cent
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.2 7 -100.72 -86.055 57.358 -114.72
## m2.1 8 -102.28 -85.525 59.140 -118.28 3.5643 1 0.05903 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m2.2)## Linear mixed model fit by REML ['lmerMod']
## Formula: wing.velocity ~ size_pc1 + Treatment.order + percLoad_cent +
## percLoad2 + (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: -87.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.89330 -0.56516 0.07513 0.49055 1.93807
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 0.004108 0.06410
## Residual 0.006188 0.07866
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.486075 0.022031 67.45
## size_pc1 -0.039588 0.009322 -4.25
## Treatment.order2 -0.008667 0.020552 -0.42
## percLoad_cent -0.380720 0.035078 -10.85
## percLoad2 -0.169751 0.138432 -1.23
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 prcLd_
## size_pc1 -0.029
## Trtmnt.rdr2 -0.381 -0.001
## percLod_cnt 0.054 0.084 0.080
## percLoad2 -0.540 0.048 -0.141 -0.150# renaming model to simplify later typing
m2 <- m2.2
##### LRTs for main effects
## size
m3 <- update(m2, .~. - Treatment.order)
anova(m2, m3, test = "Chi") # drop trt order## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m3: wing.velocity ~ size_pc1 + percLoad_cent + percLoad2 + (1 | Bee.ID)
## m2: wing.velocity ~ size_pc1 + Treatment.order + percLoad_cent +
## m2: percLoad2 + (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m3 6 -102.52 -89.954 57.260 -114.52
## m2 7 -100.72 -86.055 57.358 -114.72 0.1952 1 0.6587summary(m3)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## wing.velocity ~ size_pc1 + percLoad_cent + percLoad2 + (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: -93.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.94992 -0.54108 0.04463 0.51763 1.87287
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 0.004216 0.06493
## Residual 0.005997 0.07744
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.482722 0.020255 73.20
## size_pc1 -0.039609 0.009329 -4.25
## percLoad_cent -0.379848 0.034460 -11.02
## percLoad2 -0.179915 0.135098 -1.33
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 prcLd_
## size_pc1 -0.031
## percLod_cnt 0.090 0.083
## percLoad2 -0.643 0.048 -0.140# trt order
m4 <- update(m3, .~. - percLoad2)
anova(m3, m4) # drop percLoad## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m4: wing.velocity ~ size_pc1 + percLoad_cent + (1 | Bee.ID)
## m3: wing.velocity ~ size_pc1 + percLoad_cent + percLoad2 + (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m4 5 -102.80 -92.331 56.401 -112.80
## m3 6 -102.52 -89.954 57.260 -114.52 1.7177 1 0.19summary(m4)## Linear mixed model fit by REML ['lmerMod']
## Formula: wing.velocity ~ size_pc1 + percLoad_cent + (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: -93.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.05997 -0.53298 -0.03019 0.59699 2.09586
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 0.003737 0.06113
## Residual 0.006384 0.07990
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.465370 0.015198 96.42
## size_pc1 -0.038985 0.009136 -4.27
## percLoad_cent -0.385121 0.035065 -10.98
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1
## size_pc1 0.000
## percLod_cnt 0.000 0.095# LRT for size
m5 <- update(m4, .~. - size_pc1)
anova(m4, m5) # keep size## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m5: wing.velocity ~ percLoad_cent + (1 | Bee.ID)
## m4: wing.velocity ~ size_pc1 + percLoad_cent + (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m5 4 -89.649 -81.271 48.824 -97.649
## m4 5 -102.803 -92.331 56.401 -112.803 15.154 1 9.908e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m6 <- update(m4, .~. - percLoad_cent)
anova(m4, m6) # keep percent load## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m6: wing.velocity ~ size_pc1 + (1 | Bee.ID)
## m4: wing.velocity ~ size_pc1 + percLoad_cent + (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m6 4 -48.945 -40.567 28.472 -56.945
## m4 5 -102.803 -92.331 56.401 -112.803 55.858 1 7.791e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# summarize final model for paper
summary(m4)## Linear mixed model fit by REML ['lmerMod']
## Formula: wing.velocity ~ size_pc1 + percLoad_cent + (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: -93.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.05997 -0.53298 -0.03019 0.59699 2.09586
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 0.003737 0.06113
## Residual 0.006384 0.07990
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.465370 0.015198 96.42
## size_pc1 -0.038985 0.009136 -4.27
## percLoad_cent -0.385121 0.035065 -10.98
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1
## size_pc1 0.000
## percLod_cnt 0.000 0.095# write output
summary(m4)$coefficients ## Estimate Std. Error t value
## (Intercept) 1.46537020 0.015197520 96.421669
## size_pc1 -0.03898547 0.009135556 -4.267444
## percLoad_cent -0.38512098 0.035064635 -10.983174write.csv( summary(m4)$coefficients, file = "WingVelocity_percLoad.csv" )Wing Velocity diagnostics
# qq plot
qqnorm(resid(m4), main = "")
qqline(resid(m4)) # good
# cook's distance
plot(influence(m4, obs = TRUE), which = 'cook')
# residual plot
plot(fitted(m4), resid(m4), xlab = "fitted", ylab = "residuals")
abline(0,0)
# QQPlot for group-level effects
qqnorm(ranef(m4)$Bee.ID[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(m4)$Bee.ID[[1]]) # looks good
Wing Velocity Visualization
preds1 <- predict(m4, newdf, re.form = NA) #predict for mean of random effect
# plot of prediction for average sized bee (line)
# raw data plotted as points
ggplot(bdta, aes(x = percLoad_cent + mean(bdta$percLoad) , y = wing.velocity)) +
geom_point(alpha = 0.3) +
geom_line(data = newdf, aes(x = percLoad_cent + mean(bdta$percLoad), y = preds1)) +
labs(x = "Mass during flight (% of starved mass in g)")
Part II: Which flight measurements are most closely associated with
respiratory rate in bees – wingbeat frequency, stroke amplitude, or
wing velocity?
Multilevel model approach to see the relatedness of the response variables
- Wingbeat frequency
- Stroke amplitude
- Wing velocity
- Resp. Rate
# predicting avg resp rate, including stroke amplitude, wing velocity and wingbeat frequency
# looks like we have to remove wing velocity
vif.mer(lmer(av.resp..CO2.mL.hr. ~ strokeAmp_cent + wingVel_cent + wbf_cent + percLoad_cent + Treatment.order + size_pc1 + (1|Bee.ID), data = bdta))## strokeAmp_cent wingVel_cent wbf_cent percLoad_cent
## 190.642035 133.592345 11.644262 10.853222
## Treatment.order2 size_pc1
## 1.262757 1.364694vif.mer(lmer(av.resp..CO2.mL.hr. ~ strokeAmp_cent + wbf_cent + percLoad_cent + Treatment.order + size_pc1 + (1|Bee.ID), data = bdta)) # still worryingly high for percLoad and stroke amplitude## strokeAmp_cent wbf_cent percLoad_cent Treatment.order2
## 7.741480 1.830209 8.182520 1.255652
## size_pc1
## 1.355159mm1 <- lmer(av.resp..CO2.mL.hr. ~ (strokeAmp_cent + wbf_cent +
percLoad_cent + Treatment.order +
size_pc1)^2 +
strokeAmp2 + wbf2 + percLoad2 +
(1|Bee.ID), data = bdta)## Warning: Some predictor variables are on very different scales: consider
## rescalingsummary(mm1) # warning says to rescale## Linear mixed model fit by REML ['lmerMod']
## Formula:
## av.resp..CO2.mL.hr. ~ (strokeAmp_cent + wbf_cent + percLoad_cent +
## Treatment.order + size_pc1)^2 + strokeAmp2 + wbf2 + percLoad2 +
## (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 218.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.55796 -0.28776 -0.05692 0.20517 2.08071
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.1866 1.0893
## Residual 0.2396 0.4895
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 10.265659 0.344190 29.826
## strokeAmp_cent -0.116525 0.032532 -3.582
## wbf_cent 0.160723 0.036796 4.368
## percLoad_cent 8.410281 1.344833 6.254
## Treatment.order2 0.052521 0.163024 0.322
## size_pc1 1.936582 0.189350 10.228
## strokeAmp2 0.007777 0.002598 2.994
## wbf2 -0.002007 0.003538 -0.567
## percLoad2 12.861989 4.588437 2.803
## strokeAmp_cent:wbf_cent 0.005223 0.004034 1.295
## strokeAmp_cent:percLoad_cent -0.712136 0.201023 -3.543
## strokeAmp_cent:Treatment.order2 0.093952 0.038242 2.457
## strokeAmp_cent:size_pc1 0.020108 0.014717 1.366
## wbf_cent:percLoad_cent -0.225045 0.186731 -1.205
## wbf_cent:Treatment.order2 -0.070432 0.036113 -1.950
## wbf_cent:size_pc1 -0.042046 0.018814 -2.235
## percLoad_cent:Treatment.order2 -3.626144 1.421979 -2.550
## percLoad_cent:size_pc1 -0.467014 0.594009 -0.786
## Treatment.order2:size_pc1 -0.575909 0.150526 -3.826##
## Correlation matrix not shown by default, as p = 19 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it## fit warnings:
## Some predictor variables are on very different scales: consider rescaling# rescale numeric predictors
bdta_scaled = scale(bdta[, c("av.resp..CO2.mL.hr.",
"strokeAmp_cent",
"wbf_cent",
"percLoad_cent",
"size_pc1",
"strokeAmp2",
"wbf2",
"percLoad2")],
center = FALSE,
scale = TRUE)
colnames(bdta_scaled) <- paste0(colnames(bdta_scaled), "_scaled")
scatterplotMatrix(bdta_scaled)
bdta = cbind(bdta, bdta_scaled)
mm1 <- lmer(av.resp..CO2.mL.hr. ~ (strokeAmp_cent_scaled + wbf_cent_scaled +
percLoad_cent_scaled + Treatment.order +
size_pc1_scaled)^2 +
strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled +
(1|Bee.ID), data = bdta)
summary(mm1) # warning says to rescale## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ (strokeAmp_cent_scaled + wbf_cent_scaled +
## percLoad_cent_scaled + Treatment.order + size_pc1_scaled)^2 +
## strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 164.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.55796 -0.28776 -0.05692 0.20517 2.08071
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.1866 1.0893
## Residual 0.2396 0.4895
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 10.26566 0.34419 29.826
## strokeAmp_cent_scaled -1.43098 0.39951 -3.582
## wbf_cent_scaled 1.22494 0.28044 4.368
## percLoad_cent_scaled 2.63390 0.42117 6.254
## Treatment.order2 0.05252 0.16302 0.322
## size_pc1_scaled 3.26359 0.31910 10.228
## strokeAmp2_scaled 1.68151 0.56164 2.994
## wbf2_scaled -0.16967 0.29913 -0.567
## percLoad2_scaled 1.61846 0.57738 2.803
## strokeAmp_cent_scaled:wbf_cent_scaled 0.48889 0.37760 1.295
## strokeAmp_cent_scaled:percLoad_cent_scaled -2.73884 0.77312 -3.543
## strokeAmp_cent_scaled:Treatment.order2 1.15377 0.46962 2.457
## strokeAmp_cent_scaled:size_pc1_scaled 0.41615 0.30457 1.366
## wbf_cent_scaled:percLoad_cent_scaled -0.53715 0.44570 -1.205
## wbf_cent_scaled:Treatment.order2 -0.53679 0.27524 -1.950
## wbf_cent_scaled:size_pc1_scaled -0.54003 0.24164 -2.235
## percLoad_cent_scaled:Treatment.order2 -1.13562 0.44533 -2.550
## percLoad_cent_scaled:size_pc1_scaled -0.24648 0.31350 -0.786
## Treatment.order2:size_pc1_scaled -0.97054 0.25367 -3.826##
## Correlation matrix not shown by default, as p = 19 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it# LRT's
mm2 <- update(mm1, .~. - percLoad_cent_scaled:size_pc1_scaled)
anova(mm1, mm2)## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## mm2: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## mm2: percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## mm2: strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## mm2: Bee.ID) + strokeAmp_cent_scaled:wbf_cent_scaled + strokeAmp_cent_scaled:percLoad_cent_scaled +
## mm2: strokeAmp_cent_scaled:Treatment.order + strokeAmp_cent_scaled:size_pc1_scaled +
## mm2: wbf_cent_scaled:percLoad_cent_scaled + wbf_cent_scaled:Treatment.order +
## mm2: wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order +
## mm2: Treatment.order:size_pc1_scaled
## mm1: av.resp..CO2.mL.hr. ~ (strokeAmp_cent_scaled + wbf_cent_scaled +
## mm1: percLoad_cent_scaled + Treatment.order + size_pc1_scaled)^2 +
## mm1: strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## mm1: Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mm2 20 172.68 214.56 -66.338 132.68
## mm1 21 172.85 216.83 -65.425 130.85 1.8248 1 0.1767summary(mm2)## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## Bee.ID) + strokeAmp_cent_scaled:wbf_cent_scaled + strokeAmp_cent_scaled:percLoad_cent_scaled +
## strokeAmp_cent_scaled:Treatment.order + strokeAmp_cent_scaled:size_pc1_scaled +
## wbf_cent_scaled:percLoad_cent_scaled + wbf_cent_scaled:Treatment.order +
## wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order +
## Treatment.order:size_pc1_scaled
## Data: bdta
##
## REML criterion at convergence: 164.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.64217 -0.33793 -0.05746 0.20194 2.17217
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.1457 1.0704
## Residual 0.2449 0.4949
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 10.26035 0.34280 29.931
## strokeAmp_cent_scaled -1.38524 0.39805 -3.480
## wbf_cent_scaled 1.25452 0.27883 4.499
## percLoad_cent_scaled 2.58315 0.41822 6.177
## Treatment.order2 0.05135 0.16419 0.313
## size_pc1_scaled 3.29213 0.31545 10.436
## strokeAmp2_scaled 1.62369 0.56341 2.882
## wbf2_scaled -0.21584 0.29498 -0.732
## percLoad2_scaled 1.50825 0.56769 2.657
## strokeAmp_cent_scaled:wbf_cent_scaled 0.35059 0.33972 1.032
## strokeAmp_cent_scaled:percLoad_cent_scaled -2.57308 0.75769 -3.396
## strokeAmp_cent_scaled:Treatment.order2 1.10875 0.47154 2.351
## strokeAmp_cent_scaled:size_pc1_scaled 0.20581 0.15000 1.372
## wbf_cent_scaled:percLoad_cent_scaled -0.33822 0.37227 -0.909
## wbf_cent_scaled:Treatment.order2 -0.52811 0.27740 -1.904
## wbf_cent_scaled:size_pc1_scaled -0.57156 0.23678 -2.414
## percLoad_cent_scaled:Treatment.order2 -1.09799 0.44528 -2.466
## Treatment.order2:size_pc1_scaled -0.98374 0.25484 -3.860##
## Correlation matrix not shown by default, as p = 18 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itmm3 <- update(mm2, .~. - wbf_cent_scaled:percLoad_cent_scaled)
anova(mm2,mm3 )## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## mm3: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## mm3: percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## mm3: strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## mm3: Bee.ID) + strokeAmp_cent_scaled:wbf_cent_scaled + strokeAmp_cent_scaled:percLoad_cent_scaled +
## mm3: strokeAmp_cent_scaled:Treatment.order + strokeAmp_cent_scaled:size_pc1_scaled +
## mm3: wbf_cent_scaled:Treatment.order + wbf_cent_scaled:size_pc1_scaled +
## mm3: percLoad_cent_scaled:Treatment.order + Treatment.order:size_pc1_scaled
## mm2: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## mm2: percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## mm2: strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## mm2: Bee.ID) + strokeAmp_cent_scaled:wbf_cent_scaled + strokeAmp_cent_scaled:percLoad_cent_scaled +
## mm2: strokeAmp_cent_scaled:Treatment.order + strokeAmp_cent_scaled:size_pc1_scaled +
## mm2: wbf_cent_scaled:percLoad_cent_scaled + wbf_cent_scaled:Treatment.order +
## mm2: wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order +
## mm2: Treatment.order:size_pc1_scaled
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mm3 19 172.39 212.18 -67.193 134.39
## mm2 20 172.68 214.56 -66.338 132.68 1.7113 1 0.1908summary(mm3)## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## Bee.ID) + strokeAmp_cent_scaled:wbf_cent_scaled + strokeAmp_cent_scaled:percLoad_cent_scaled +
## strokeAmp_cent_scaled:Treatment.order + strokeAmp_cent_scaled:size_pc1_scaled +
## wbf_cent_scaled:Treatment.order + wbf_cent_scaled:size_pc1_scaled +
## percLoad_cent_scaled:Treatment.order + Treatment.order:size_pc1_scaled
## Data: bdta
##
## REML criterion at convergence: 165
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.64387 -0.33501 -0.04398 0.30860 2.16600
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.1202 1.0584
## Residual 0.2491 0.4991
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 10.32638 0.33546 30.783
## strokeAmp_cent_scaled -1.24881 0.37193 -3.358
## wbf_cent_scaled 1.24684 0.27979 4.456
## percLoad_cent_scaled 2.45007 0.39328 6.230
## Treatment.order2 0.08464 0.16095 0.526
## size_pc1_scaled 3.22835 0.30701 10.515
## strokeAmp2_scaled 1.37824 0.50568 2.726
## wbf2_scaled -0.32636 0.27002 -1.209
## percLoad2_scaled 1.14468 0.41345 2.769
## strokeAmp_cent_scaled:wbf_cent_scaled 0.06994 0.14601 0.479
## strokeAmp_cent_scaled:percLoad_cent_scaled -2.12997 0.59831 -3.560
## strokeAmp_cent_scaled:Treatment.order2 1.00405 0.46268 2.170
## strokeAmp_cent_scaled:size_pc1_scaled 0.24015 0.14555 1.650
## wbf_cent_scaled:Treatment.order2 -0.52241 0.27921 -1.871
## wbf_cent_scaled:size_pc1_scaled -0.56075 0.23620 -2.374
## percLoad_cent_scaled:Treatment.order2 -0.96560 0.42207 -2.288
## Treatment.order2:size_pc1_scaled -0.92574 0.24880 -3.721##
## Correlation matrix not shown by default, as p = 17 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itmm4 <- update(mm3, .~. - strokeAmp_cent_scaled:size_pc1_scaled)
anova(mm3, mm4) # keep interaction## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## mm4: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## mm4: percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## mm4: strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## mm4: Bee.ID) + strokeAmp_cent_scaled:wbf_cent_scaled + strokeAmp_cent_scaled:percLoad_cent_scaled +
## mm4: strokeAmp_cent_scaled:Treatment.order + wbf_cent_scaled:Treatment.order +
## mm4: wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order +
## mm4: Treatment.order:size_pc1_scaled
## mm3: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## mm3: percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## mm3: strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## mm3: Bee.ID) + strokeAmp_cent_scaled:wbf_cent_scaled + strokeAmp_cent_scaled:percLoad_cent_scaled +
## mm3: strokeAmp_cent_scaled:Treatment.order + strokeAmp_cent_scaled:size_pc1_scaled +
## mm3: wbf_cent_scaled:Treatment.order + wbf_cent_scaled:size_pc1_scaled +
## mm3: percLoad_cent_scaled:Treatment.order + Treatment.order:size_pc1_scaled
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mm4 18 175.30 213.00 -69.651 139.30
## mm3 19 172.39 212.18 -67.193 134.39 4.916 1 0.02661 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(mm3)## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## Bee.ID) + strokeAmp_cent_scaled:wbf_cent_scaled + strokeAmp_cent_scaled:percLoad_cent_scaled +
## strokeAmp_cent_scaled:Treatment.order + strokeAmp_cent_scaled:size_pc1_scaled +
## wbf_cent_scaled:Treatment.order + wbf_cent_scaled:size_pc1_scaled +
## percLoad_cent_scaled:Treatment.order + Treatment.order:size_pc1_scaled
## Data: bdta
##
## REML criterion at convergence: 165
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.64387 -0.33501 -0.04398 0.30860 2.16600
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.1202 1.0584
## Residual 0.2491 0.4991
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 10.32638 0.33546 30.783
## strokeAmp_cent_scaled -1.24881 0.37193 -3.358
## wbf_cent_scaled 1.24684 0.27979 4.456
## percLoad_cent_scaled 2.45007 0.39328 6.230
## Treatment.order2 0.08464 0.16095 0.526
## size_pc1_scaled 3.22835 0.30701 10.515
## strokeAmp2_scaled 1.37824 0.50568 2.726
## wbf2_scaled -0.32636 0.27002 -1.209
## percLoad2_scaled 1.14468 0.41345 2.769
## strokeAmp_cent_scaled:wbf_cent_scaled 0.06994 0.14601 0.479
## strokeAmp_cent_scaled:percLoad_cent_scaled -2.12997 0.59831 -3.560
## strokeAmp_cent_scaled:Treatment.order2 1.00405 0.46268 2.170
## strokeAmp_cent_scaled:size_pc1_scaled 0.24015 0.14555 1.650
## wbf_cent_scaled:Treatment.order2 -0.52241 0.27921 -1.871
## wbf_cent_scaled:size_pc1_scaled -0.56075 0.23620 -2.374
## percLoad_cent_scaled:Treatment.order2 -0.96560 0.42207 -2.288
## Treatment.order2:size_pc1_scaled -0.92574 0.24880 -3.721##
## Correlation matrix not shown by default, as p = 17 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itmm5 <- update(mm3, .~. - wbf_cent_scaled:Treatment.order)
anova(mm3, mm5) # keep interaction## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## mm5: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## mm5: percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## mm5: strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## mm5: Bee.ID) + strokeAmp_cent_scaled:wbf_cent_scaled + strokeAmp_cent_scaled:percLoad_cent_scaled +
## mm5: strokeAmp_cent_scaled:Treatment.order + strokeAmp_cent_scaled:size_pc1_scaled +
## mm5: wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order +
## mm5: Treatment.order:size_pc1_scaled
## mm3: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## mm3: percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## mm3: strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## mm3: Bee.ID) + strokeAmp_cent_scaled:wbf_cent_scaled + strokeAmp_cent_scaled:percLoad_cent_scaled +
## mm3: strokeAmp_cent_scaled:Treatment.order + strokeAmp_cent_scaled:size_pc1_scaled +
## mm3: wbf_cent_scaled:Treatment.order + wbf_cent_scaled:size_pc1_scaled +
## mm3: percLoad_cent_scaled:Treatment.order + Treatment.order:size_pc1_scaled
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mm5 18 176.32 214.02 -70.159 140.32
## mm3 19 172.39 212.18 -67.193 134.39 5.9309 1 0.01488 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(mm3)## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## Bee.ID) + strokeAmp_cent_scaled:wbf_cent_scaled + strokeAmp_cent_scaled:percLoad_cent_scaled +
## strokeAmp_cent_scaled:Treatment.order + strokeAmp_cent_scaled:size_pc1_scaled +
## wbf_cent_scaled:Treatment.order + wbf_cent_scaled:size_pc1_scaled +
## percLoad_cent_scaled:Treatment.order + Treatment.order:size_pc1_scaled
## Data: bdta
##
## REML criterion at convergence: 165
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.64387 -0.33501 -0.04398 0.30860 2.16600
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.1202 1.0584
## Residual 0.2491 0.4991
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 10.32638 0.33546 30.783
## strokeAmp_cent_scaled -1.24881 0.37193 -3.358
## wbf_cent_scaled 1.24684 0.27979 4.456
## percLoad_cent_scaled 2.45007 0.39328 6.230
## Treatment.order2 0.08464 0.16095 0.526
## size_pc1_scaled 3.22835 0.30701 10.515
## strokeAmp2_scaled 1.37824 0.50568 2.726
## wbf2_scaled -0.32636 0.27002 -1.209
## percLoad2_scaled 1.14468 0.41345 2.769
## strokeAmp_cent_scaled:wbf_cent_scaled 0.06994 0.14601 0.479
## strokeAmp_cent_scaled:percLoad_cent_scaled -2.12997 0.59831 -3.560
## strokeAmp_cent_scaled:Treatment.order2 1.00405 0.46268 2.170
## strokeAmp_cent_scaled:size_pc1_scaled 0.24015 0.14555 1.650
## wbf_cent_scaled:Treatment.order2 -0.52241 0.27921 -1.871
## wbf_cent_scaled:size_pc1_scaled -0.56075 0.23620 -2.374
## percLoad_cent_scaled:Treatment.order2 -0.96560 0.42207 -2.288
## Treatment.order2:size_pc1_scaled -0.92574 0.24880 -3.721##
## Correlation matrix not shown by default, as p = 17 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itmm6 <- update(mm3, .~. - strokeAmp_cent_scaled:wbf_cent_scaled)
anova(mm3, mm6) # drop interaction## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## mm6: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## mm6: percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## mm6: strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## mm6: Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order +
## mm6: strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order +
## mm6: wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order +
## mm6: Treatment.order:size_pc1_scaled
## mm3: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## mm3: percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## mm3: strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## mm3: Bee.ID) + strokeAmp_cent_scaled:wbf_cent_scaled + strokeAmp_cent_scaled:percLoad_cent_scaled +
## mm3: strokeAmp_cent_scaled:Treatment.order + strokeAmp_cent_scaled:size_pc1_scaled +
## mm3: wbf_cent_scaled:Treatment.order + wbf_cent_scaled:size_pc1_scaled +
## mm3: percLoad_cent_scaled:Treatment.order + Treatment.order:size_pc1_scaled
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mm6 18 171.04 208.74 -67.518 135.04
## mm3 19 172.39 212.18 -67.193 134.39 0.6498 1 0.4202summary(mm6)## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order +
## strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order +
## wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order +
## Treatment.order:size_pc1_scaled
## Data: bdta
##
## REML criterion at convergence: 163.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6717 -0.3034 -0.0680 0.3528 2.2000
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.1046 1.0510
## Residual 0.2438 0.4937
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 10.33167 0.33224 31.097
## strokeAmp_cent_scaled -1.21537 0.36143 -3.363
## wbf_cent_scaled 1.20177 0.26092 4.606
## percLoad_cent_scaled 2.45005 0.38932 6.293
## Treatment.order2 0.06881 0.15588 0.441
## size_pc1_scaled 3.18673 0.29199 10.914
## strokeAmp2_scaled 1.36399 0.49926 2.732
## wbf2_scaled -0.27572 0.24625 -1.120
## percLoad2_scaled 1.12004 0.40566 2.761
## strokeAmp_cent_scaled:percLoad_cent_scaled -2.08682 0.58439 -3.571
## strokeAmp_cent_scaled:Treatment.order2 0.96049 0.44821 2.143
## strokeAmp_cent_scaled:size_pc1_scaled 0.19023 0.10033 1.896
## wbf_cent_scaled:Treatment.order2 -0.48672 0.26616 -1.829
## wbf_cent_scaled:size_pc1_scaled -0.51630 0.21509 -2.400
## percLoad_cent_scaled:Treatment.order2 -0.98064 0.41651 -2.354
## Treatment.order2:size_pc1_scaled -0.88736 0.23299 -3.809##
## Correlation matrix not shown by default, as p = 16 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itmm7 <- update(mm6, .~. - wbf_cent_scaled:size_pc1_scaled)
anova(mm6, mm7) # keep## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## mm7: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## mm7: percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## mm7: strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## mm7: Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order +
## mm7: strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order +
## mm7: percLoad_cent_scaled:Treatment.order + Treatment.order:size_pc1_scaled
## mm6: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## mm6: percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## mm6: strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## mm6: Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order +
## mm6: strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order +
## mm6: wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order +
## mm6: Treatment.order:size_pc1_scaled
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mm7 17 176.71 212.32 -71.357 142.71
## mm6 18 171.04 208.74 -67.518 135.04 7.6766 1 0.005594 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(mm6)## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order +
## strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order +
## wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order +
## Treatment.order:size_pc1_scaled
## Data: bdta
##
## REML criterion at convergence: 163.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6717 -0.3034 -0.0680 0.3528 2.2000
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.1046 1.0510
## Residual 0.2438 0.4937
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 10.33167 0.33224 31.097
## strokeAmp_cent_scaled -1.21537 0.36143 -3.363
## wbf_cent_scaled 1.20177 0.26092 4.606
## percLoad_cent_scaled 2.45005 0.38932 6.293
## Treatment.order2 0.06881 0.15588 0.441
## size_pc1_scaled 3.18673 0.29199 10.914
## strokeAmp2_scaled 1.36399 0.49926 2.732
## wbf2_scaled -0.27572 0.24625 -1.120
## percLoad2_scaled 1.12004 0.40566 2.761
## strokeAmp_cent_scaled:percLoad_cent_scaled -2.08682 0.58439 -3.571
## strokeAmp_cent_scaled:Treatment.order2 0.96049 0.44821 2.143
## strokeAmp_cent_scaled:size_pc1_scaled 0.19023 0.10033 1.896
## wbf_cent_scaled:Treatment.order2 -0.48672 0.26616 -1.829
## wbf_cent_scaled:size_pc1_scaled -0.51630 0.21509 -2.400
## percLoad_cent_scaled:Treatment.order2 -0.98064 0.41651 -2.354
## Treatment.order2:size_pc1_scaled -0.88736 0.23299 -3.809##
## Correlation matrix not shown by default, as p = 16 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itmm8 <- update(mm6, .~. - percLoad_cent_scaled:Treatment.order)
anova(mm6, mm8) # keep## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## mm8: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## mm8: percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## mm8: strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## mm8: Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order +
## mm8: strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order +
## mm8: wbf_cent_scaled:size_pc1_scaled + Treatment.order:size_pc1_scaled
## mm6: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## mm6: percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## mm6: strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## mm6: Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order +
## mm6: strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order +
## mm6: wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order +
## mm6: Treatment.order:size_pc1_scaled
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mm8 17 176.86 212.47 -71.431 142.86
## mm6 18 171.04 208.74 -67.518 135.04 7.8259 1 0.00515 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(mm6)## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order +
## strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order +
## wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order +
## Treatment.order:size_pc1_scaled
## Data: bdta
##
## REML criterion at convergence: 163.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6717 -0.3034 -0.0680 0.3528 2.2000
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.1046 1.0510
## Residual 0.2438 0.4937
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 10.33167 0.33224 31.097
## strokeAmp_cent_scaled -1.21537 0.36143 -3.363
## wbf_cent_scaled 1.20177 0.26092 4.606
## percLoad_cent_scaled 2.45005 0.38932 6.293
## Treatment.order2 0.06881 0.15588 0.441
## size_pc1_scaled 3.18673 0.29199 10.914
## strokeAmp2_scaled 1.36399 0.49926 2.732
## wbf2_scaled -0.27572 0.24625 -1.120
## percLoad2_scaled 1.12004 0.40566 2.761
## strokeAmp_cent_scaled:percLoad_cent_scaled -2.08682 0.58439 -3.571
## strokeAmp_cent_scaled:Treatment.order2 0.96049 0.44821 2.143
## strokeAmp_cent_scaled:size_pc1_scaled 0.19023 0.10033 1.896
## wbf_cent_scaled:Treatment.order2 -0.48672 0.26616 -1.829
## wbf_cent_scaled:size_pc1_scaled -0.51630 0.21509 -2.400
## percLoad_cent_scaled:Treatment.order2 -0.98064 0.41651 -2.354
## Treatment.order2:size_pc1_scaled -0.88736 0.23299 -3.809##
## Correlation matrix not shown by default, as p = 16 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itmm9 <- update(mm6, .~. - Treatment.order:size_pc1_scaled)
anova(mm6, mm9) # keep## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## mm9: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## mm9: percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## mm9: strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## mm9: Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order +
## mm9: strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order +
## mm9: wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order
## mm6: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## mm6: percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## mm6: strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## mm6: Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order +
## mm6: strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order +
## mm6: wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order +
## mm6: Treatment.order:size_pc1_scaled
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mm9 17 186.96 222.57 -76.481 152.96
## mm6 18 171.04 208.74 -67.518 135.04 17.926 1 2.296e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(mm6)## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order +
## strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order +
## wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order +
## Treatment.order:size_pc1_scaled
## Data: bdta
##
## REML criterion at convergence: 163.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6717 -0.3034 -0.0680 0.3528 2.2000
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.1046 1.0510
## Residual 0.2438 0.4937
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 10.33167 0.33224 31.097
## strokeAmp_cent_scaled -1.21537 0.36143 -3.363
## wbf_cent_scaled 1.20177 0.26092 4.606
## percLoad_cent_scaled 2.45005 0.38932 6.293
## Treatment.order2 0.06881 0.15588 0.441
## size_pc1_scaled 3.18673 0.29199 10.914
## strokeAmp2_scaled 1.36399 0.49926 2.732
## wbf2_scaled -0.27572 0.24625 -1.120
## percLoad2_scaled 1.12004 0.40566 2.761
## strokeAmp_cent_scaled:percLoad_cent_scaled -2.08682 0.58439 -3.571
## strokeAmp_cent_scaled:Treatment.order2 0.96049 0.44821 2.143
## strokeAmp_cent_scaled:size_pc1_scaled 0.19023 0.10033 1.896
## wbf_cent_scaled:Treatment.order2 -0.48672 0.26616 -1.829
## wbf_cent_scaled:size_pc1_scaled -0.51630 0.21509 -2.400
## percLoad_cent_scaled:Treatment.order2 -0.98064 0.41651 -2.354
## Treatment.order2:size_pc1_scaled -0.88736 0.23299 -3.809##
## Correlation matrix not shown by default, as p = 16 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itmm10 <- update(mm6 , .~. - strokeAmp_cent_scaled:Treatment.order)
anova(mm6, mm10) # keep interaction## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## mm10: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## mm10: percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## mm10: strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## mm10: Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:size_pc1_scaled +
## mm10: wbf_cent_scaled:Treatment.order + wbf_cent_scaled:size_pc1_scaled +
## mm10: percLoad_cent_scaled:Treatment.order + Treatment.order:size_pc1_scaled
## mm6: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## mm6: percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## mm6: strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## mm6: Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order +
## mm6: strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order +
## mm6: wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order +
## mm6: Treatment.order:size_pc1_scaled
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mm10 17 176.19 211.79 -71.093 142.19
## mm6 18 171.04 208.74 -67.518 135.04 7.1495 1 0.007499 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mm11 <- update(mm6, .~. - strokeAmp_cent_scaled:Treatment.order)
anova(mm6, mm11) # keep interaction## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## mm11: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## mm11: percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## mm11: strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## mm11: Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:size_pc1_scaled +
## mm11: wbf_cent_scaled:Treatment.order + wbf_cent_scaled:size_pc1_scaled +
## mm11: percLoad_cent_scaled:Treatment.order + Treatment.order:size_pc1_scaled
## mm6: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## mm6: percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## mm6: strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## mm6: Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order +
## mm6: strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order +
## mm6: wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order +
## mm6: Treatment.order:size_pc1_scaled
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mm11 17 176.19 211.79 -71.093 142.19
## mm6 18 171.04 208.74 -67.518 135.04 7.1495 1 0.007499 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(mm6)## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order +
## strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order +
## wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order +
## Treatment.order:size_pc1_scaled
## Data: bdta
##
## REML criterion at convergence: 163.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6717 -0.3034 -0.0680 0.3528 2.2000
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.1046 1.0510
## Residual 0.2438 0.4937
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 10.33167 0.33224 31.097
## strokeAmp_cent_scaled -1.21537 0.36143 -3.363
## wbf_cent_scaled 1.20177 0.26092 4.606
## percLoad_cent_scaled 2.45005 0.38932 6.293
## Treatment.order2 0.06881 0.15588 0.441
## size_pc1_scaled 3.18673 0.29199 10.914
## strokeAmp2_scaled 1.36399 0.49926 2.732
## wbf2_scaled -0.27572 0.24625 -1.120
## percLoad2_scaled 1.12004 0.40566 2.761
## strokeAmp_cent_scaled:percLoad_cent_scaled -2.08682 0.58439 -3.571
## strokeAmp_cent_scaled:Treatment.order2 0.96049 0.44821 2.143
## strokeAmp_cent_scaled:size_pc1_scaled 0.19023 0.10033 1.896
## wbf_cent_scaled:Treatment.order2 -0.48672 0.26616 -1.829
## wbf_cent_scaled:size_pc1_scaled -0.51630 0.21509 -2.400
## percLoad_cent_scaled:Treatment.order2 -0.98064 0.41651 -2.354
## Treatment.order2:size_pc1_scaled -0.88736 0.23299 -3.809##
## Correlation matrix not shown by default, as p = 16 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it# move on to squared terms
mm12 <- update(mm6 , .~. - wbf2_scaled)
anova(mm6, mm11) # keep wbf2## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## mm11: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## mm11: percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## mm11: strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## mm11: Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:size_pc1_scaled +
## mm11: wbf_cent_scaled:Treatment.order + wbf_cent_scaled:size_pc1_scaled +
## mm11: percLoad_cent_scaled:Treatment.order + Treatment.order:size_pc1_scaled
## mm6: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## mm6: percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## mm6: strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## mm6: Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order +
## mm6: strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order +
## mm6: wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order +
## mm6: Treatment.order:size_pc1_scaled
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mm11 17 176.19 211.79 -71.093 142.19
## mm6 18 171.04 208.74 -67.518 135.04 7.1495 1 0.007499 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(mm6)## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order +
## strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order +
## wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order +
## Treatment.order:size_pc1_scaled
## Data: bdta
##
## REML criterion at convergence: 163.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6717 -0.3034 -0.0680 0.3528 2.2000
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.1046 1.0510
## Residual 0.2438 0.4937
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 10.33167 0.33224 31.097
## strokeAmp_cent_scaled -1.21537 0.36143 -3.363
## wbf_cent_scaled 1.20177 0.26092 4.606
## percLoad_cent_scaled 2.45005 0.38932 6.293
## Treatment.order2 0.06881 0.15588 0.441
## size_pc1_scaled 3.18673 0.29199 10.914
## strokeAmp2_scaled 1.36399 0.49926 2.732
## wbf2_scaled -0.27572 0.24625 -1.120
## percLoad2_scaled 1.12004 0.40566 2.761
## strokeAmp_cent_scaled:percLoad_cent_scaled -2.08682 0.58439 -3.571
## strokeAmp_cent_scaled:Treatment.order2 0.96049 0.44821 2.143
## strokeAmp_cent_scaled:size_pc1_scaled 0.19023 0.10033 1.896
## wbf_cent_scaled:Treatment.order2 -0.48672 0.26616 -1.829
## wbf_cent_scaled:size_pc1_scaled -0.51630 0.21509 -2.400
## percLoad_cent_scaled:Treatment.order2 -0.98064 0.41651 -2.354
## Treatment.order2:size_pc1_scaled -0.88736 0.23299 -3.809##
## Correlation matrix not shown by default, as p = 16 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itmm13 <- update(mm6, .~. - percLoad2_scaled)
anova(mm6, mm13) # keep## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## mm13: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## mm13: percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## mm13: strokeAmp2_scaled + wbf2_scaled + (1 | Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled +
## mm13: strokeAmp_cent_scaled:Treatment.order + strokeAmp_cent_scaled:size_pc1_scaled +
## mm13: wbf_cent_scaled:Treatment.order + wbf_cent_scaled:size_pc1_scaled +
## mm13: percLoad_cent_scaled:Treatment.order + Treatment.order:size_pc1_scaled
## mm6: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## mm6: percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## mm6: strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## mm6: Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order +
## mm6: strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order +
## mm6: wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order +
## mm6: Treatment.order:size_pc1_scaled
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mm13 17 179.69 215.29 -72.844 145.69
## mm6 18 171.04 208.74 -67.518 135.04 10.652 1 0.001099 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mm14 <- update(mm6, .~. - strokeAmp2_scaled)
anova(mm6, mm14) # keep## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## mm14: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## mm14: percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## mm14: wbf2_scaled + percLoad2_scaled + (1 | Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled +
## mm14: strokeAmp_cent_scaled:Treatment.order + strokeAmp_cent_scaled:size_pc1_scaled +
## mm14: wbf_cent_scaled:Treatment.order + wbf_cent_scaled:size_pc1_scaled +
## mm14: percLoad_cent_scaled:Treatment.order + Treatment.order:size_pc1_scaled
## mm6: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## mm6: percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## mm6: strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## mm6: Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order +
## mm6: strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order +
## mm6: wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order +
## mm6: Treatment.order:size_pc1_scaled
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mm14 17 179.11 214.71 -72.555 145.11
## mm6 18 171.04 208.74 -67.518 135.04 10.073 1 0.001504 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(mm6) # final model for paper## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ strokeAmp_cent_scaled + wbf_cent_scaled +
## percLoad_cent_scaled + Treatment.order + size_pc1_scaled +
## strokeAmp2_scaled + wbf2_scaled + percLoad2_scaled + (1 |
## Bee.ID) + strokeAmp_cent_scaled:percLoad_cent_scaled + strokeAmp_cent_scaled:Treatment.order +
## strokeAmp_cent_scaled:size_pc1_scaled + wbf_cent_scaled:Treatment.order +
## wbf_cent_scaled:size_pc1_scaled + percLoad_cent_scaled:Treatment.order +
## Treatment.order:size_pc1_scaled
## Data: bdta
##
## REML criterion at convergence: 163.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6717 -0.3034 -0.0680 0.3528 2.2000
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.1046 1.0510
## Residual 0.2438 0.4937
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 10.33167 0.33224 31.097
## strokeAmp_cent_scaled -1.21537 0.36143 -3.363
## wbf_cent_scaled 1.20177 0.26092 4.606
## percLoad_cent_scaled 2.45005 0.38932 6.293
## Treatment.order2 0.06881 0.15588 0.441
## size_pc1_scaled 3.18673 0.29199 10.914
## strokeAmp2_scaled 1.36399 0.49926 2.732
## wbf2_scaled -0.27572 0.24625 -1.120
## percLoad2_scaled 1.12004 0.40566 2.761
## strokeAmp_cent_scaled:percLoad_cent_scaled -2.08682 0.58439 -3.571
## strokeAmp_cent_scaled:Treatment.order2 0.96049 0.44821 2.143
## strokeAmp_cent_scaled:size_pc1_scaled 0.19023 0.10033 1.896
## wbf_cent_scaled:Treatment.order2 -0.48672 0.26616 -1.829
## wbf_cent_scaled:size_pc1_scaled -0.51630 0.21509 -2.400
## percLoad_cent_scaled:Treatment.order2 -0.98064 0.41651 -2.354
## Treatment.order2:size_pc1_scaled -0.88736 0.23299 -3.809##
## Correlation matrix not shown by default, as p = 16 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it# write output
summary(mm9)$coefficients ## Estimate Std. Error
## (Intercept) 10.50939334 0.3719988
## strokeAmp_cent_scaled -0.78693814 0.4224585
## wbf_cent_scaled 0.79953361 0.2846480
## percLoad_cent_scaled 1.99184379 0.4507705
## Treatment.order2 0.03758213 0.2003348
## size_pc1_scaled 2.62560994 0.2567266
## strokeAmp2_scaled 0.88098584 0.6199180
## wbf2_scaled -0.39713261 0.3031080
## percLoad2_scaled 0.74256129 0.5006662
## strokeAmp_cent_scaled:percLoad_cent_scaled -1.30791510 0.7164999
## strokeAmp_cent_scaled:Treatment.order2 0.69012764 0.5682134
## strokeAmp_cent_scaled:size_pc1_scaled 0.17030338 0.1302434
## wbf_cent_scaled:Treatment.order2 0.18640327 0.2572330
## wbf_cent_scaled:size_pc1_scaled -0.45047772 0.2380927
## percLoad_cent_scaled:Treatment.order2 -0.86679594 0.5231384
## t value
## (Intercept) 28.2511470
## strokeAmp_cent_scaled -1.8627583
## wbf_cent_scaled 2.8088498
## percLoad_cent_scaled 4.4187536
## Treatment.order2 0.1875966
## size_pc1_scaled 10.2272624
## strokeAmp2_scaled 1.4211327
## wbf2_scaled -1.3102017
## percLoad2_scaled 1.4831465
## strokeAmp_cent_scaled:percLoad_cent_scaled -1.8254227
## strokeAmp_cent_scaled:Treatment.order2 1.2145571
## strokeAmp_cent_scaled:size_pc1_scaled 1.3075783
## wbf_cent_scaled:Treatment.order2 0.7246475
## wbf_cent_scaled:size_pc1_scaled -1.8920268
## percLoad_cent_scaled:Treatment.order2 -1.6569152write.csv( summary(mm11)$coefficients, file = "resp_oth_Coefs_percLoad.csv" )met rate with other predictors diagnostics
# qq plot
qqnorm(resid(mm6), main = "")
qqline(resid(mm6)) # good, but two outliers
# residual plot
plot(fitted(mm6), resid(mm6), xlab = "fitted", ylab = "residuals")
abline(0,0)
# QQPlot for group-level effects
qqnorm(ranef(mm6)$Bee.ID[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(mm6)$Bee.ID[[1]]) # looks good
# visualize model:
sjp.lmer(mm6, type = 'fe')## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.
met rate with other variable visualization
newdf <- data.frame(strokeAmp_cent_scaled = 0,
wbf_cent_scaled = 0,
percLoad_cent_scaled = seq(min(bdta$percLoad_cent_scaled),
max(bdta$percLoad_cent_scaled),
length.out = 100),
Treatment.order = factor(1),
size_pc1_scaled = 0,
strokeAmp2_scaled = 0,
wbf2_scaled = 0
)
newdf$percLoad2_scaled <- scale(newdf$percLoad_cent_scaled^2)
preds1 <- predict(mm6, newdf, re.form = NA)
# plot of prediction for average sized bee (line)
# raw data plotted as points
ggplot(bdta, aes(x = percLoad_cent_scaled , y = av.resp..CO2.mL.hr.)) +
geom_point(alpha = 0.3) +
geom_line(data = newdf, aes(x = percLoad_cent_scaled, y = preds1)) +
labs(x = "Mass during flight (% of starved mass in g)")
# visualize metabolic rate vs. wbf
newdf <- data.frame(strokeAmp_cent_scaled = 0,
wbf_cent_scaled = seq(min(bdta$wbf_cent_scaled),
max(bdta$wbf_cent_scaled),
length.out = 100),
percLoad_cent_scaled = 0,
percLoad2_scaled = 0,
Treatment.order = factor(1),
size_pc1_scaled = 0,
strokeAmp2_scaled = 0
)
newdf$wbf2_scaled <- scale(newdf$wbf_cent_scaled^2)
preds1 <- predict(mm6, newdf, re.form = NA)
# plot of prediction for average sized bee (line)
# raw data plotted as points
ggplot(bdta, aes(x = wbf_cent_scaled , y = av.resp..CO2.mL.hr.)) +
geom_point(alpha = 0.3) +
geom_line(data = newdf, aes(x = wbf_cent_scaled, y = preds1)) +
labs(x = "Wing beat frequency (centered and scaled)")
# visualize metabolic rate vs. stroke
newdf <- data.frame(
wbf_cent_scaled = 0,
wbf2_scaled = 0,
percLoad_cent_scaled = 0,
percLoad2_scaled = 0,
Treatment.order = factor(1),
size_pc1_scaled = 0,
strokeAmp_cent_scaled = seq(min(bdta$strokeAmp_cent_scaled),
max(bdta$strokeAmp_cent_scaled),
length.out = 100)
)
newdf$strokeAmp2_scaled <- scale(newdf$strokeAmp_cent_scaled^2)
preds1 <- predict(mm6, newdf, re.form = NA)
# plot of prediction for average sized bee (line)
# raw data plotted as points
ggplot(bdta, aes(x = strokeAmp_cent_scaled , y = av.resp..CO2.mL.hr.)) +
geom_point(alpha = 0.3) +
geom_line(data = newdf, aes(x = strokeAmp_cent_scaled, y = preds1)) +
labs(x = "Stroke Amplitude (centered and scaled)")
This model is fundamentally different from the previous one. It’s saying, when we account for differences in wingbeat freq, stroke amplitude, and bee size, we sill see an effect of percent load on metabolic rate.
OR, we could say that holding all other variables constant, an increase in wingbeat frequency is associated with an increase in metabolic rate.